Chaotic behavior and damage spreading in the Glauber Ising model - a master 

equation approach 
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We investigate the sensitivity of the time evolution of a kinetic Ising model with Glauber dynamics 
against the initial conditions. To do so we apply the "damage spreading" method, i.e., we study 
the simultaneous evolution of two identical systems subjected to the same thermal noise. We derive 
a master equation for the joint probability distribution of the two systems. We then solve this 
master equation within an effective-field approximation which goes beyond the usual mean-field 
approximation by retaining the fluctuations though in a quite simplistic manner. The resulting 
effective-field theory is then applied to different physical situations. It is used to analyze the fixed 
points of the master equation and their stability and to identify regular and chaotic phases of the 
Glauber Ising model. We also discuss the relation of our results to directed percolation. 



PACS numbers: 05.40. +j, 64.60. Ht, 75.40. Gb 
I. INTRODUCTION 

The physics of dynamic phase transitions and dynamic 
critical phenomena has been a subject of great interest 
for the last two decades. Whereas the dynamic behavior 
at and close to usual static phase transitions is well un- 
derstood much less is known about dynamic phase 
transitions which do not have a static counterpart. Some- 
times it is not even known whether or not a particular 
dynamic transition coincides with an equilibrium phase 
transition. 

One of these dynamic phenomena is the so-called 
"damage spreading" The central question of this 

problem is how a small perturbation (called the damage) 
in a cooperative system changes during the further time 
evolution. Among the simplest of such cooperative sys- 
tems are kinetic Ising models where the above question 
has been investigated by means of Monte-Carlo simula- 
tions |^,^| . In these simulations two identical Ising mod- 
els with different initial conditions are subjected to the 
same thermal noise, i.e. the same random numbers are 
used in the Monte-Carlo procedure. In analogy to the 
physics of chaotic dynamics [|| the differences in the mi- 
croscopic configurations of the two systems are then used 
to characterize the dynamics and to distinguish regular 
and chaotic phases, depending on external parameters 
(e.g. temperature, magnetic field). 

Later the name "damage spreading" has also been ap- 
plied to a different though related type of investigations 
in which the two systems are not identical. Instead, one 
or several spins in one of the copies are permanently fixed 
in one direction. Therefore the equilibrium properties of 
the two systems are different and the microscopic dif- 
ferences between the two copies can be related to static 
and dynamic correlation functions . Note that in this 
type of simulations it is not essential to use identical noise 
(i.e. random numbers) for the two systems. Instead it is 
only a convenient method to reduce the statistical error. 

Whereas this second type of damage spreading is well 



understood and established as a method to numerically 
calculate equilibrium properties, much less is known 
about the original problem of damage spreading, viz. 
how sensitive is the dynamics of the Ising model to differ- 
ent initial conditions. In particular, there are no rigorous 
results on the transition between regular and chaotic be- 
havior (called the "spreading transition"). 

There are two different mechanisms by which the dam- 
age can spread in a kinetic Ising model. First, the dam- 
age can spread within a single ergodic component (i.e. a 
pure state or free energy valley) of the system. This is 
the case for Glauber or Metropolis dynamics. Numerical 
simulations here consistently give a transition tempera- 
ture slightly lower than the equilibrium critical tempera- 
ture 0. Grassberger conjectured that the spreading 
transition falls into the universality class of directed per- 
colation if it does not coincide with another phase tran- 
sition. This was supported by high-precision numerical 
simulations for the Glauber Ising model pd[ . 

Second, the damage can spread when the system se- 
lects one of the free energy valleys at random after a 
quench from high temperatures to below the equilibrium 
critical temperature. This is the only mechanism to pro- 
duce damage spreading in an Ising model with heat-bath 
dynamics. In this case the spreading temperature seems 
to coincide with the equilibrium critical temperature be- 
low which the two pure states separate [f4 |12|]l3|l . Thus, 
at the spreading point there are long-range static corre- 
lations in the systems, and the transition is expected to 
fall into a universality class different from directed per- 
colation. 

In this paper we investigate the damage spreading in 
the Glauber Ising model by deriving and solving a mas- 
ter equation for the time evolution of a joint probability 
distribution for two identical systems with different ini- 
tial conditions and subjected to the same thermal noise. 
The paper is organized as follows. In Section IIA we 
define the model. Transition probabilities between the 
states of a spin pair are calculated in section IIB and 
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the master equation for the joint probability distribution 
is derived in section IIC. We discuss how to construct 
a mean-field approximation for this equation in section 
IIIA. In sections IIIB, IIIC and HID we present solutions 
of the master equation within this approximation for dif- 
ferent physical situations. Finally, section IV is dedicated 
to conclusions and an outlook on future work. A short 
account of part of this work has already been published 
JLlj together with a comparison to the heat-bath Ising 
model. 



II. A MASTER EQUATION FOR DAMAGE 
SPREADING 

A. The Glauber Ising model 

We consider two identical kinetic Ising models with N 
sites described by the Hamiltonians and given 
by 



H 



(n) 



■> X! ■ / '-- s 



(n) g(ri) 



in) 



(1) 



where is an Ising variable with the values ±1 and 
n = 1, 2 distinguishes the two copies. Jy is the exchange 
interaction between the spins and h denotes an external 
magnetic field. The dynamics of the Ising models is given 
by the Glauber algorithm, i.e. in every time step a lattice 
site i is chosen at random (the same site in both copies) . 
The new value of the spin at this site is calculated ac- 
cording to 



S^(t + 1) = sign\ v[hl n '(t)]-- + Sr(t) 



1 



(2) 



where the transition probality v(x) is given by the usual 
Glauber expression 



v(x) 



e x/T /(e */T + g-x/T) 



(3) 



Here (t) = J2j J*j Sj (*) + h is the local magnetic 
field at site i and (discretized) time t in the system n. 
^i(t) £ [0, 1) is a random number which is identical for 
both systems, and T denotes the temperature. The spins 
at all sites different from site i are unchanged within this 
time step. 

The central quantity in any damage spreading process 
is the distance between the two systems in phase space, 
called the Hamming distance (or the damage). It is de- 
fined by 



D(t) 



^IX ) (t)-fiP , (t)i 



2N 



(4) 



and identical to the portion of sites where the spins in 
the two systems differ. 



In order to describe the simultaneous time evolution of 
the two systems H^> and we define a variable v{t) 
at each lattice site which describes the state of a spin 
pair (S^\S^). It has the values v = ++ for = 
= 1, +- for S« = = 1, -+ for -S« = 

= 1 and for = = -1. A complete 

configuration of the two Ising models is thus described 
by the set . . . , vn}- 

Since we are interested in the time evolution not for a 
single sequence of £i(t), but in ^-averaged quantities we 
consider a whole ensemble of system pairs (H l - 1 \H < - 2 * > ) 
and define a probability distribution 

P(v 1 ,...,v N ,t) = /j2H 5 « i Mt)) ( 5 ) 

\v t (t) i I 

where (•) denotes the ensemble average. 

B. Transition probabilities 

In order to formulate a master equation for the prob- 
ability distribution P(vi, . . . , vn , t) we need to know the 
transition probabilities w(y — > [/,) between the states v of 
a spin pair. Since the Glauber dynamics (2) involves only 
a single lattice site within each time step, we have to con- 
sider transitions between the states v of a single site only. 

Let us look, e.g., at the transition of site i from state 

to ++. This corresponds to both and chang- 
ing from —1 to 1. According to the Glauber dynamics 
(2) this requires v(h. 



& > and v{hf ] ) - & > 0. 
Since v(h) is a monotonous function of h both equations 

are simultaneously fulfilled for u[min.(/ij , h^)] — £j > 0. 
Because £j is a random number taken from a uniform 
distribution between and 1, the transition probability 
is given by 



v[mux(h (1 \h^)]. 



Analogously, for a transition from state 



two inequalities v(h 



0-h 



£i > and v(h. 



(6) 

to H — the 

6 < o 

have to be fulfilled. Since v(h) is a monotonous function 

of h this is only possible for > hf^ . The transition 
probability is obviously given by 



w{- 



+-) = 6(fr (1) - M 2) ) \v{h^) - v(h,W)] . (7) 



The transition probabilities w(v — -> /i) fulfill the fol- 
lowing symmetry relations 



w(+-\ > v) — w( > v) 

w(-\ *■!/)= w( 1 ► v) 



(8a) 
(8b) 



for any state v as can easily be seen by making the sub- 



stitutions S$ n \t) 



-S^ n \t) and &(t) 



on 



the right-hand side of (2). 

The remaining transition probabilities can be calcu- 
lated along the same lines as above. They are summa- 
rized in Table I. 
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TABLE I. Transition probabilities w{y — > n) between the 
states of a spin pair 









v[— max 










++ 


v[min(h ' , h )] 









-+ 




'v{h {2) ) -v(h w ) 








+- 


Q(h^-h^) 


[v(h (1) ) - v(h i2) ) 




-+ 






e(-h.w -h^)[v(~h^)-v(h^)] 


-+ 




++ 


e(hw + h^) \v(h^) - v(-h,W)] 


-+ 




-+ 


w[min(- 






-+ 




+- 


v[— max( 







C. The master equation 

Having calculated the transition probabilities between 
the states v of a spin pair we are now in the position 
to write down the equation of motion for the probalitity 
distribution P(i>\, ...fjvi t). It has the form of a usual 
master equation 

^P(v x ,...,v N ,t) = 

N 
N 

+ Y P{vu---,Vi,---,v N ,t)w(ni^>Vi) (9) 

i— 1 fii^Ui 

where the term in the second line describes the de- 
crease of P(vi, ...v N ,t) due to the initial configuration 
{v\, . . . vm} being changed at one of the sites i from Vi to 
Hi. The term in the third line of the master equation de- 
scribes the increase of P(vi, ...v N ,t) due to "scattering" 
from all the other states into {v\, . . . vn}. Note that we 
have suppressed the factor 1/N in the transition proba- 
bilities which corresponds to random selection of one of 
the lattice sites in every time step. This neglect corre- 
sponds to a redefinition of the time scale (which is now 
independent of the system size) and does not change the 
dynamic behavior. 

This master equation contains, of course, the full dif- 
ficulty of the dynamic many-body problem. A complete 
solution is therefore out of question, and one has to resort 
to approximation methods. In the following section we 
discuss how to construct a mean-field like approximation 
to the master equation (9). 

III. EFFECTIVE-FIELD APPROXIMATION 

Usually a mean-field theory of a phase transition can 
be obtained by taking the range of the interaction to 
infinity: 

Jij = J /N for all i,j (10) 



In the thermodynamic limit N — > oo this suppresses all 

(n) 

fluctuations. In particular, the local magnetic fields h\ 
of all sites in one system become equal and identical to 
the mean-field value J^m. Since the two Ising models 
pW and H (2) are thermodynamically identical this leads 
to h^p — h[ 2 K However, some of the transition proba- 
bilities depend on the existence of fluctuations (see table 
I), i.e. w{v — > fi) g° to zero with M 1 ) — — ► 0. In 

particular, this is true for w( ► — h) and w( ► 

H — ) which are responsible for increasing the damage D. 
Consequently, if the thermodynamic limit and the limit 
of infinite range of the interaction are taken at a too 
early stage of the calculation, the resulting model does 
not show any spreading of the damage. To circumvent 
these problems we develop a slightly more sophisticated 
effective-field approximation that retains the fluctuations 
though in a quite simplistic manner. As will be shown in 
section IIIC, taking the range of the interaction to infin- 
ity within the framework of this approximation yields a 
sensible limit. 



A. Effective-field theory for short-range models 

The central idea of this effective-field method is to 
retain the fluctuations but to treat the fluctuations 
at different sites as statistically independent. This 
amounts to approximating the probability distribution 
P(vi, . . . ,i/N,t) by a product of identical single-site dis- 
tributions P v , 

N 

P(v u ...,v N ,t) = l[P Vi (t)- (11) 

1=1 

Inserting this into the master equation (9) and summing 
over all states of sites i = 2 ... N gives an equation of 
motion for the single-site distribution P Vl , 

jP Vl = J2 l-PxiW("i - Mi) + P„ 1 W(hi - i/i)], 

(12) 

where 

W(fii -» t^.) = (w(fxi -» vi)) P (13) 

is the transition probability averaged over the states 
of all sites i ^ 1 according to the distribution P Vi . Since 
all sites of the systems are equivalent the site index i will 
be suppressed from now on. 

Note that the average magnetizations m^, of the 
two systems and the Hamming distance D can be easily 
expressed in terms of P„, 

m« = P ++ +P+_ -P_+ -P__, (14a) 
m (2) = P++ - P+_ + P_+ - P__, (14b) 
P = P+-+P-+. (14c) 
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TABLE II. Probabilities for the states of the three neigh- 
boring sites and resulting local magnetic fields and h 1 ^ 1 
for the two-dimensional case on the hexagonal lattice 



h {2) 


3J 


J 


-J 


-3J 


3J 


p3 


3PI+P-+ 


3P++P^+ 


P 3 - + 


J 


3PZ+P+- 


3PI+P— + 


3P+_P£+ + 


3P- + P!_ 






6P++P+-P-+ 


6P+-P-+P-- 




-J 


3P++PI- 


3P|_P_+ + 


3P++P£_ + 


3P-+Pl_ 






6P++P+-P— 


6P+-P-+P-- 




-3J 


Pl- 


3P+_P— 


3P+_P^_ 


P 3 .. 



So far the considerations have been rather general, in 
the following subsections we will apply the general for- 
malism to different physical situations. In section IIIB 
we investigate a two-dimensional system with short-range 
interactions and vanishing external field. We determine 
not only the location of the spreading transition but also 
calculate the stationary states of the systems. Section 
IIIC deals with the limit of infinite-range interactions, 
and in sec. HID we study the influence of an external 
magnetic field on the spreading transition. 

B. Solution of a two-dimensional model 

In this subsection we investigate the damage spreading 
for a two-dimensional Glauber Ising model on a hexag- 
onal lattice (with each site having three nearest neigh- 
bors). The interaction is taken to be a nearest-neighbor 
interaction of strength J and the external magnetic field 
is set to zero. 

In order to solve the master equation (12) for the 
single-site distribution P v we first determine the effec- 
tive transition probabilities W(y — > fi). Let us calculate 
the probabilities for a particular site i. To this end we 
have to average the transition probabilities given in ta- 
ble I with respect to the states of all other sites of the 
system. However, since the interaction is between near- 
est neighbors, the transition probabilities depend on the 
states of these three neighbors of site i only. Each of 
the neighbors can be in one of four states, thus we have 
to consider 64 different configurations of the neighboring 
sites. The probabilities for these configurations and the 
resulting local magnetic fields are given in table II. 

With the help of tables I and II the averaged transition 
probabilities (13) can be easily calculated by adding up 
the contributions of all 64 configurations. The resulting 
expression are quite lengthy though simple. Therefore 
we present only the example 

W( ► ++) = 

(w( ► ++)) = («[min(fcW,fc( 2) ]) = 

Pl + v 3 + 3Pl + P- +Vl + 3P++J^ + «_i + Pl + v- 3 + 



(3P2 + p + _ + 3Pl + p__ + QP ++ P + _P_ + ) Vl + 
(3P+_p2 + + 6P + _p_ + p__) v _ 1 + 3P_ + p2_u_ 3 + 
(3P ++ P£_ + 3P^_P_ + + 6P ++ P + _P__)v_ 1 + 
(3P++Pl_ + 6P+_P_+P__>_i + 3P-+Pl_v- 3 + 
{Pl_ + 3P|_P__ + 3P+_Pf_ + P!_)w- 3 (15) 

with 

v n = v(nJ). (16) 

Equations of motion for the magnetizations and 
rrS 2 ^ as well as for the damage D can be derived by insert- 
ing the definitions (14) into the single-site master equa- 
tion (12). After some manipulations the equations of 
motion for the magnetizations read 

j t m {n) = m (n) j-1 + | [tanh(3 J/T) + tanh( J/T)]\ 

+ \{m {n) f [tanh(3 J/T) - 3 tanh( J/T)] . (17) 

These equations are, of course, identical to the equation 
of motion of the magnetization derived for a single system 
within the same framework of statistically independent 
fluctuations. The point at which the coefficient of the 
term linear in m on the right-hand side of (17) changes 
sign defines the (equilibrium) critical temperature Tq of 
the Ising model within our approximation. Tq is thus 
determined by 

^[tanh(3J/T c ) + tanh( J/T c )] = 1 (18) 

which gives Tq/J ~ 2.104. The stationary solution of 
(17) can be used to determine the magnetization as a 
function of temperature. For temperatures T < Tc we 
obtain 

{m M ? I Nnh(3 J/T) + tanh( J/T)} 1 

1 ; | tanh( J/T) - \ tanh(3 J/T) ' 1 ' 

We now turn to the discussion of the Hamming dis- 
tance D. After inserting (14) into (12) the equation of 
motion of the Hamming distance D can be written as 

jD = (1 - D)[W( ► +-) + W( ► -+)] + 

+ D[-l + W(-+^+-) + W(-+^-+)}. (20) 

Since in the following we will be mainly interested in 
the stationary solutions of this equation we restrict the 
considerations to cases where both systems are in equi- 
librium at the beginning of the damage spreading pro- 
cess. In doing so we exclude, however, all phenomena 
connected with the behavior after a quench from high 
temperatures to temperatures below Tc- These phenom- 
ena require an investigation of the early time behavior 
and will be analyzed elsewhere |H| . 
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It is now useful to distinguish three cases, (i) damage 
spreading in the paramagnetic phase (T > Tc), (ii) the 
ferromagnetic phase (T < Tc) where both systems are 
in the same pure state (i.e., free energy valley), mA 1 ' — 
m/ 2 ^ — m and (iii) the ferromagnetic phase (T < Tc) 
where the two systems are in different pure states, m^ 1 ' — 
—vn/ 2 " 1 = m. 



1. Paramagnetic phase 

In the paramagnetic phase all P v can be expressed in 
terms of D: 



P ++ =P- = \{l-D), 

Pa- = P-+ = Id. 



(21a) 
(21b) 



By inserting this into the transition probabilities W{y —> 
Ii) calculated from (13) and table II the equation of mo- 
tion (20) of the Hamming distance D can be written as 



d \, 
—D = -ID 

dt 2 y 



3D 2 + 2D 3 ) tanh(3 J/T). 



(22) 



This equation has three stationary solutions (fixed 
points), viz. D\ = which corresponds to both sys- 
tems being identical, D| = 1 where S*' 1 ' = — for all 
sites and D% = 1/2 which corresponds to the two sys- 
tems being completely uncorrelated |lq ]. To determine 
the stability of the fixed points we linearize the equation 
of motion (22) in dk = D — D^. The linearized equation 
has the solution 



d k (t) = rffeo 



„A fc t 



(23) 



with Ai = A 2 = \ tanh(3 J/T) and A 3 = -\ tanh(3 J/T). 
Consequently, the only stable fixed point is Dg = 1/2. In 
the whole paramagnetic phase the damage spreads and 
asymptotically reaches the value D — 1/2. If the two 
systems start very close together (D small initially) their 
distance in phase space increases exponentially with a 
Lyapunov exponent Ai = ^tanh(3J/T). Therefore the 
Glauber dynamics shows chaotic behavior in the whole 
paramagnetic phase. Note, that for large temperatures 
the Lyapunov exponent Ai goes to zero as Ai ~ 3J/T. 
Thus, the time it takes the system to reach the station- 
ary state D = D3 = 1/2 diverges for T — > oo. This 
has recently also been found in simulations jl7| . The de- 
pendence of the Lyapunov exponent on temperature is 
presented in fig. 1. 



2. Ferromagnetic phase with rw 1 ' = = m 

In this paragraph we study the case where both sys- 
tems are in the same free energy valley. The single-site 
probabilities P v can be expressed in terms of D and m: 
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FIG. 1. Magnetization m, asymptotic Hamming distance 
D* and Lyapunov exponent Ai as functions of temperature 
for the Glauber Ising model. Below T c the curve for D* has 
two branches corresponding to the two systems being in the 
same or in different free energy valleys. 



P++ = i(l-D + m), 
P— = \{l-D-m), 



P+- = P-+ = ~D. 

+ + 2 



(24a) 
(24b) 
(24c) 



After inserting this into the averaged transition probabil- 
ities (13) the equation of motion of the Hamming distance 
takes the form 

^-D = i(D - 3D 2 + 2D 3 ) tanh(3J/T) 

- ^TO 2 [2L>tanh( J/T) 

-D 2 tanh( J/T) + D 2 tanh(3 J/T)]. (25) 

This equation has two fixed points D* in the interval 
[0,1]. The first fixed point is D\ — 0. By linearizing 
(25) in d-y = D — D\ we investigate the stability of this 
fixed point. We again find that d\{t) follows the exponen- 
tial law (23) with Ai = itanh(3 J/T) -\m 2 tanh( J/T). 
Using for m 2 the expression (19) it is easy to discuss the 
behavior of Ai. For temperatures larger than a spreading 
temperature T$ which is defined by 



3m 2 tanh( J/T s ) = tanh(3 J/T s ) 



(26) 



the Lyapunov exponent Ai is positive and thus the fixed 
point D1 is unstable. For T < T$ the Lyapunov ex- 
ponent Ai is negative and the fixed point D\ is stable. 
Consequently, the Glauber dynamics is chaotic for tem- 
peratures above T$ but regular below. Eq. (26) gives 
T s ~ 1.739 J w 0.826T C . 

For temperatures T > T$ the equation of motion (25) 
possesses another fixed point Dg with < D^ < 1/2 
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which is always stable. Its temperature dependence is 
presented in fig. 1. Close to the spreading temperature 
the asymptotic Hamming distance D3 increases linearly 
with T — T s which corresponds to the spreading transi- 
tion being of 2nd order. The order parameter exponent 
(3, defined by D* = |T-T S |' 3 is given by (3 =. In contrast 
to the paramagnetic phase, where the two systems even- 
tually become completely uncorrelated, for T s < T < T c 
the asymptotic Hamming distance D is always smaller 
than 1/2 so that the two systems remain partially cor- 
related (as it must be the case since both systems are 
in the same free energy valley). Directly at the spread- 
ing point the term linear in D in (25) vanishes. For 
small Hamming distances the equation of motion now 
reads dD/dt oc — D 2 which gives a power-law behavior 
D(t) oc t~ s . The critical exponent is given by 5 = 1. 



3. Ferromagnetic phase with m 



(i) 



,(2) 



We now turn to the case where the two systems are in 
different free energy valleys. The single-site probabilities 
P v can be expressed in terms of D and m: 



P. 



++ 



p- = la-D), 



p+- 



-(D + m), 



P_ + = -(D-m). 



(27a) 
(27b) 
(27c) 



With this substitutions the equation of motion (20) of 
the Hamming distance can be written as 



1 



dt 



D = -{D — 3D + 2D A ) tanh(3 J/T) 



2 
3 

+ -to 2 [tanh(3 J/T) + tanh( J/T) 



+D 2 tanh(3 J/T) 



2D tanh(3 J/T) 
D 2 tanh( J/T)}. (28) 



Analogously to the preceeding paragraph, this equation 
possesses two fixed points. The fixed point D\ = 1 exits 
for all temperatures. It is stable for temperatures below 
Ts and unstable above. For T > Ts (28) has another 
fixed point , D\ with 1/2 < D\ < 1 which is always 
stable. Its temperature dependence is given in fig. 1. 



the limit z — > 00. To obtain a physically sensible limit 
we scale the interaction strength with z, J = Jq/z. 

In the limit z — > 00 the thermodynamics is described 
by the usual mean-field theory. The equilibrium critical 
temperature is given by Tc = Jo and in the ferromagnetic 
phase the magnetization is determined by the equation 
of state 



to = tanh(mJ /T)- 



(29) 



In order to determine the spreading temperature Ts it 
is sufficient to study the equation of motion (20) of the 
Hamming distance to linear order in D. To this end we 

have to determine W(— H > H — ) and W( — I > — h) to 

zeroth order in D but W{ ► H — ) and W( ► — h) 

to linear order in D. 

To zeroth order in D we have = = h and in 
the limit z — > 00 h is ^-distributed at h = J^m. The 
transition probabilities are thus given by (see table I) 

W{— H ► H — ) = v[— max(— h^,h^)] 

= v(-J \m\), (30a) 

W(— H > — +) = w[min(— h^,h^)] 

= v(-J \m\). (30b) 

We now calculate W{ ► H — ) and W{ ► — h) 

to linear order in D. These transition probabilities do 
not have a zeroth-order contribution. In linear order in 
D only those configurations of the z neighboring sites 
contribute for which the two systems differ in the state 
of a single site. In this case and differ by 2J /z. 
Therefore we obtain 

W{ ► +-) = 6(fr (1) - hW)[v(h^) - v(h^)} 

= z P + _(2J /z) v'(J m). (31a) 

Here v'(h) is the derivative of v with respect to its argu- 
ment. The additional factor z in the second line arises 
since each of the z neighbors can be the one where the 
two systems differ. Analogously we obtain 

W{ ► +-) = e(/t (2) - hW)[v{hW) - v{h^)} 

= z P_+(2J /z) v'(J m). (31b) 

By inserting these results for the transition probabilities 
into the equation of motion (20) of the Hamming distance 
we find 



C. The limit of high dimensions 

In this subsection we study damage spreading in the 
Glauber Ising model in the limit of high dimensions, i.e. 
in the mean-field limit proper. Within the framework of 
our effective field approach high dimensions correspond 
to high coordination numbers, i.e. high numbers of near- 
est neighbors. We therefore consider a Glauber Ising 
model on a lattice with z nearest neighbors and study 
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D = XD 



(32) 



where the Lyapunov exponent A is given by 

e^p(~J \m\/T) 



A = -1 + 2 
+ 



exp(J H/T) + exp(-JoM/T) 
4J /T 



[exp(J TO/T) + cxp(- J to/T)] 2 ' 
This can be simplified to 



(33) 
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FIG. 2. Magnetization m and Lyapunov exponent Ai as 
functions of temperature for the Glauber Ising model with 
vanishing external field in the limit of high dimensions. 



A 



(1 



)Jo/T. 



(34) 



In the paramagnetic phase (m = 0) the Lyapunov expo- 
nent is simply X = J /T > 0. Thus the Glauber dynam- 
ics is chaotic in the whole paramagnetic phase. 

The temperature dependence of the Lyapunov expo- 
nent A in the ferromagnetic phase is presented in fig. 
2. A changes sign at T s ~ 0.827J = 0.827T C . Conse- 
quently, the dynamics is chaotic for temperatures larger 
than Ts = 0.827J and regular for temperatures smaller 
than Ts- Note that the value for T$/Tc for the two- 
dimensional model of sec. IIIB is very close to but not 
identical to the value for the case z — > oo. 



■< 0.0 
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FIG. 3. Lyapunov exponent A as a function of temperature 
for external fields h — 0, 0.2, 0.4 and 0.6 (from up to down). 



nal field shifts the spreading temperature to higher val- 
ues, thus suppressing chaotic behavior and stabilizing the 
regular phase. The phase boundary between the chaotic 
and the regular phase can be easily determined by solving 
the equation A = 0. The resulting phase diagram is pre- 
sented in fig. 4. For comparison we also give simulation 
results for a three-dimensional Glauber Ising model. 
An investigation of (35) and (36) for large temperatures 
shows that the spreading temperature Ts(h) diverges for 
h/Jo^l as 

T s (h)/J = 1/(1 - h/J ) for h/J -> 1. (37) 

Consequently, for external fields h > Jq the dynamics is 
always regular. 



D. Damage spreading in a field 



IV. CONCLUSIONS 



In this subsection we generalize the effective-field the- 
ory to a finite external magnetic field h. For simplicity, 
we do this only for the model introduced in the proceed- 
ing subsection, viz. the limiting case of high dimensions. 

The equation of state (29) has to be replaced by 



m = tanh[(mJo + h)/T]. 



(35) 



Analogously, in all transition probabilities W(v — ► fi) the 
term Jqiti has to be replaced by Jqiti + h. After inserting 
the transition probabilities into the equation of motion 
(20) of the Hamming distance one finds d/dt D = XD 
and the Lyapunov exponent can again be expressed in 
terms of the magnetization: 



A = -TO+ (1 -m 2 )J a /T. 



(36) 



The temperature and field dependence of the Lyapunov 
exponent is illustrated in fig. 3. Obviously, the exter- 



To summarize, we have developed a master equation 
approach to damage spreading and applied it to the 
Glauber Ising model. The master equation is an exact 
description of the damage spreading problem, it does not 
contain any approximations. We have then solved the 
master equation within an effective-field theory for vari- 
ous physical situations. 

In this final section we discuss some aspects which have 
not been covered yet. First , we compare the results 
of our effective-field theory with numerical simulations 
of damage spreading of the Glauber Ising model in two 
and three dimensions |) 11 In agreement with the 
simulation results we find a spreading transition below 
the equilibrium critical temperature of the Ising model. 
Our mean-field value T s /T c w 0.827 is considerably lower 
than the latest numerical values of 0.992 for a two- 
dimensional and 0.922 for a three-dimensional Glauber 
Ising model. We expect our value to be exact, how- 
ever, for an infinite-dimensional model or, equivalently, 
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FIG. 4. Phase diagram of damage spreading in the Glauber 
Ising model in the limit of infinite dimensions. The full line 
shows the result of our theory, the squares are simulation 
results of Le Caer with T and h rescaled by Tc- 



An obvious idea is to include quenched disorder into 
the Hamiltonian of the Ising model either in the form of 
a random external field or in the form of random interac- 
tions. Such systems have been numerically investigated 
in some detail, in particular in the case of random inter- 
actions Recently some interesting results have also 
been achieved for random fields |l7j ]. Some investiga- 
tions on the application of the master equation approach 
to disordered systems are in progress. 
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for infinite range of the interaction. Grassberger |lC|] 
conjectured that the damage spreading transition in the 
Glauber Ising model is in the universality class of directed 
percolation. Our results are compatible with that, since 
the values of the critical exponents (3 (which describes 
the dependence of the stationary damage on the reduced 
temperature) and 8 (which describes the time decay of 
the damage at the spreading temperature) are identical 
to the mean-field values (3 = 8 = 1 of directed percola- 
tion. 

Second, we want to clarify the relation to damage 
spreading in an Ising model with heat-bath dynamics. As 
already discussed in the introduction the heat-bath Ising 
model does not show any spreading of damage within a 
single pure state (free energy valley) . When applying our 
effective-field theory to the heat-bath Ising model we find 
|L4j only a single fixed point D"{ — if both systems are 
in the same pure state M . If the two systems are in dif- 
ferent pure states (for T < Tc) we also find a single fixed 
point only, viz D\ = m. Thus there is no chaotic be- 
havior within one pure state. However, the damage can 
spread (or, at least, will not heal) in the heat-bath Ising 
model if the two copies start in different pure states or 
choose to develop into different pure states after a quench 
from high temperatures. For this case a mean-field the- 
ory similar to our's has been considered before [ p0| . 

Finally, we discuss possible extensions of the present 
theory. In principle, the master equation approach of 
sec. II can be applied to any damage spreading problem 
in which the dynamics of a single system is given by a 
stochastic map as in (2) (or a more general map that 
involves several sites in each time step). It would be 
very interesting to remap the master equation onto a field 
theory and then apply renormalization group methods to 
determine the critical behavior. 
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